Self-consistent theory of the thermal softening and instability of simple crystals 
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We consider the thermal softening of crystals due to anharmonicity. Self-consistent methods 
find a maximum temperature for a stable crystal, which gives an upper bound to the melting 
temperature. Previous workers have shown that the self-consistent harmonic approximation (SCHA) 
gives misleading results for the thermal stability of crystals. The reason is that the most important 
diagrams in the perturbation expansion around harmonic theory are not included in the SCHA. An 
alternative approach is to solve a self-consistent Dyson equation (SCA) for a selection of diagrams, 
the simplest being a (3-|-4)-SCA. However, this gives an unsatisfactory comparison to numerical 
results on the thermal and quantum melting of two-dimensional (2D) Coulomb-interacting particles 
(equivalent to vortex-lattice melting in two and three dimensions). We derive an improved self- 
consistent method, the two-vertex-SCHA, which gives much better agreement to the simulations. 
Our method allows for accurate calculation of the thermal softening of the shear modulus for 2D 
crystals and for the lattice of vortex-lines in type-H superconductors. 

PACS numbers: 63.70.+h, 64.70. Dv, 62.20.Dc, 74.60.Ec, 74.60. Ge 



I. INTRODUCTION 

It was realized long ago on theoretical grapnds that 
most melting transitions should be first order,! although 
exceptions include the possibility of |defect-mediated 
melting in two-dimensional (2D) crystalta and the mean- 
field transition between the superconducting vortexJat- 
tice and the normal state in type-II superconductors.B As 
with most first-order transitions, there is no good sim- 
ple theory of melting to compare with our theoretical 
understanding of many continuous phase transitions: A 
first-order transition occurs when the free energy of one 
phase as a function of temperature crosses the free en- 
ergy of another phase (see Fig. |l|). These two phases are 
qualitatively different and generally demand theoretical 
treatments with distinct approximations, such that a de- 
tailed comparison of free energies is not possible. The 
solid-liquid transition is a prime example of this prob- 
lem. We can treat the solid phase well within elasticity 
theory, but this is not a useful description of the liquid. 

A common approach to this problem is to concenixate 
on the solid phase and calculate its stability limito (a 
complementary method is-|to treat the liquid, e.g. within 
density functional theoryjj and find the lowest tempera- 
ture for the liquid phase to exist). A plausible mechanism 
for the crystal instability comes from anharmonic effects 
that may soften the lattice when thermal fluctuations are 
large. To this end one can devjelop a self-consistent har- 
monic approximation (SCHA) where an effective har- 
monic theory is used that better approximates the true 
anharmonic system for finite thermal fluctuations. We 
can then identify the limit to a self-consistent solution 
with the instability temperature of the solid. While it is 
unclear whether or not the instability point T„ is close 



to the melting temperature T^, it does give an upper 
bound (see Section III), and represents the first step to 
understanding the possible mechanisms of melting for a 
given system. 
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FIG. 1. Schematic diagram of the free energy as a func- 
tion of temperature for a first-order phase transition. The 
thermodynamic phase is the one of lowest free energy, and 
the possibility of a (meta-)stable solution of a phase exists 
beyond the temperature of the transition. 



An example of this hunt for the anharmonic instabil- 
ity point occurred several years ago in the context of 
2D melting. An early SCHA calculation by Platzman 
and FukuyamaEl for a 2D electron system found a dra- 
matic stiffening of the Wigner crystal as a function of 
temperature, and the anharmonic instability point was 
reached well above the dislocation-mediated melting tem- 
perature. Later, Fisher considered anharmonic correc- 
tions from a perturbation expansion around the harmonic 
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theory, and could show that the lowest-order tempera^ 
ture correction to the shear modulus was downwards.u 
The opposite result of the SCHA calculation comes from 
the fact that it misses contributions from odd anhaij^ 
monicities. A series of papers by Lozovik and coworkeraa 
emphasized this fact, and developed an alternative self- 
consistent procedure from the Dyson equation for the two 
lowest order diagrams of the expansion. This alternative 
method was named the (3-|-4)-self consistent approxima- 
tion (SCA). (The SCA is distinguished from the SCHA 
in that self-consistency comes from the Dyson resumma- 
tion of a set of skeletal diagrams in the perturbation the- 
ory, rather than by using an effective harmonic theory to 
take averages over the true Hamiltonian.) An important 
result of Lozovik's work was that for long-range interac- 
tions in 2D, the anharmonic instability point lies below 
the dislocation mediated melting temperature, implying 
a first-order rather than continuous transition. 

Our own interest in the melting problem originates in 
studies of the vortex system in high-Tc superconductors 
(see Ref. |l^ for a review.) There is excellent experimental 
evidence for a first-order melting transition of the vortps 
lattice in clean crystals of RiaStpCaCuzOg (BSCCO)tll 
and YBa2Cu307_5 (YBCO),Eltj with a latent heat and 
magnetization step consistent with our theoretical mod- 
els of the vortex system.O However, the detailed under- 
standing of where and why this transition j-takes place 
has so far relied on numerical simulations .E3~l3 With this 
in mind, we have tried to apply the above analytic ex^ 
tensions to elasticity theory within the London modelEiJ 
to derive the instability point of the vortex crystal. By 
comparing our results to the numerical simulations of 
Rcf. |l| we have found, in agreement with Lozovik, that 
the SCHA seriously underestimat es the thermal softening 
of the vortex lattice (see Section VIC). However, appli- 
cation of the diagrammatic (3+4)-SCA did not lead to a 
good quantitative comparison, and so we have developed 
an improved self-consistent method that includes all di- 
agrams in the (3-|-4)-SCA and the SCHA. Because the 
improved method uses an effective harmonic theory, but 
is equivalent to the Dyson equation for all two- vertex di- 
agrams, we call it the two-vertex-SCHA. The successful 
comparison of this new method to numerical results for 
the vortex lattice forms the main result of this paper. 

The article is structured as follows. We first define the 
harmonic approximation for a generalized crystal. Then 
in Section III we describe the self-consistent harmonic 



approximation and its justification within a variational 
approach. In Section IV we review the perturbation ex- 
pansion about the harmonic approximation, look at the 
diagrams corresponding to the SCHA and see why this 
approximation can give such misleading results. The al- 
ternative (3-|-4)-SCA of Lozovik et al. is also introduced, 
and we discuss the advantages and disadvantages over 
the SCHA. After learning the diagrammatic details of 
these self-consistent approximations, in Section ^ we in- 
troduce our improved self-consistent method, the two- 
vertex-SCHA, that includes all diagrams in both the 



SCHA, (3-|-4)-SCA and more, while keeping some of the 
simplicity in treatment of the SCHA. Finally in Sec- 
tion VI we apply this improved method to some phys- 
First we consider the thermal melting in 
"one-component plasma" (with 



ical problems, 
two dimensions of the 
ln(i?) interactions) and the Wigner crystal (with 1/R 
interactions). In contrast to the (3-|-4)-SCA, our re- 
sults from the two-vertex-SCHA leave open the possi- 
bility of dislocation-mediated melting. We also calculate 
the quantum melting of 2D bosons, equivalent to ther- 
mal melting of vortex lines. We find a significant thermal 
softening of the vortex lattice- which could be responsible 
for the observed peak effectE3 in critical current close to 
the melting transition in YBCO. 



II. REVIEW OF HARMONIC THEORY 

We consider a classical d-dimensional crystal at finite 
temperature, where the particle positions make small 
fiuctuations {u^} about their equilibrium sites {R;^}. All 
thermodynamic properties are controlled by the partition 
function. 
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(1) 



where the Hamiltonian -ff [u^] gives the energy increase 
from the ground state for a given configuration of 
displacements.E3 We work in Fourier space (which diag- 
onalizes a harmonic Hamiltonian), and define. 



(2) 



where n is the particle density. The inverse transform is 
an integral over the Brillouin zone. 
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In the harmonic approximation (HA) we expand the 
Hamiltonian to second order in displacements. 



with 
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The elastic matrix is defined by 



(4) 

(5) 
(6) 

(7) 
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Here L'' is the total volume, which determines the dis- 
cretization of fc-space. We define the harmonic Green's 
function as the inverse of the elastic matrix, GQ^Ck) = 

The HA, being a quadratic theory, is completely solv- 
able, with the size of fluctuations given by the equipar- 
tition theorem. The harmonic propagator is the thermal 
average, 

= r^'^(ki+k2)G^''(ki), (8) 

where Zq is the partition function for the harmonic 
Hamiltonian H2. For dimensions d > 2 a useful mea- 
sure of the fluctuations is the mean-square displacement 
of each particle (in 2D we need to be careful about di- 
vergences at long wavelengths, see Appendix^. This is 
just the integral over k of the propagator, 

{u\ = (^."(R)«"(R)>o =T [ ^Gr(k). (9) 

Jbz [^^) 



III. THE SELF-CONSISTENT HARMONIC 
APPROXIMATION 

From Eq. (^) we see how the mean displacements in- 
crease with temperature. At some stage the anharmonic 
corrections to Eq. (^) become important. One way to 
treat these anharmonic effects is to take the quadratic 
form as a trial Hamiltonian, 

with the elastic matrix as a set of variationaLjparame- 
ters. We can then use the general inequalitynJ for the 
free energy F = —TlnZ, 



F<Ft 



(H - Hf 



(11) 



so that for a given trial Hamiltonian we get the best ap- 
proximation to the free energy by minimizing the right 
hand side. It is straightforward to show that for the 
quadratic form this minimization is satisfied by the ef- 
fective elastic matrix. 



(12) 



where the average is over the distribution defined by the 
trial Hamiltonian ( |l0|) at the given temperature. This 
must be solved self-consistently, which defines the SCHA. 

Let us consider the general case of a system of particles 
with pairwise interactions. 



(13) 



where R^ are the ideal crystal positions with nearest 
neighbor separation a. Eq. (p2|) then becomes. 



n ^^(1 — cosk • R^) 



r=R„+u„ 



(14) 



This average over the second derivative of the potential is 
conveniently written in terms of the displacement fluctu- 
ations by Fourier transforming the interaction potential. 



(15) 
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Because the fields in a harmonic potential will have Gaus- 
sian fluctuations, we have used the property (e*'^'^) = 
g-ife-(a; ) £pj. ^ Gaussian distributed variable x. Al- 
though not necessary for calculations, the formalism sim- 
plifies greatly if we can make the approximation that 
((m^ - Uq)(u'1 - u^)) = {2/d)Saf3{u'^), which means ig- 
noring correlations between the fluctuations of different 
particle positions. In this case we flnd, 

$f (k) = ^(Q^ + k)"(Q^ -I- k)^/(|Q,, + k|) 



(16) 



where are the reciprocal lattice vectors, and f{q) = 
V{q)e~'^ {" )/d^ The fluctuations (w^) are calculated from 
the equipartition theorem (H) but with the Green's func- 
tion given as the inverse of the effective elastic matrix. 
Thus within this homogeneous approximation, we have 
turned the problem of finding (k) for all k to a self- 
consistent equation for a single parameter, (m^). From 
( p6| ) it is clear that at large {u^} the elastic matrix will 
be suppressed from the zero temperature limit, which in 
turn increases (u^) above the harmonic expectation. It 
is this feature which drives an anharmonic instability of 
the crystal phase at a certain temperature. 

The SCHA in the "independent fluctuations" approxi- 
mation has an interesting feature for the case of Coulomb 
interactions {l/r in three dimensions, Inr in two dimen- 
sions etc.). By definition the Laplacian of a Coulomb 
interaction is zero, and one can show that the renor- 
malized interaction, and therefore the effective elastic 
matrix, is a non-perturbative function of (u^). This 
is shown explicitly for the case of logarithmic interac- 
tions in two dimensions in Appendix ^. For the 2D 
case, the single particle fiuctuations diverge (see Ap- 
pendix ^) even at low temperatures where the shear 
modulus remains non-zeroB However, the shear modulus 



3 



fit = linife^^o [^t^{ky)/K'^] is dominated by the nearest- 
neighbor fluctuations (u^)ni 



i(|u(R^)-u(R^+i)|2) 



the formula for /it(('u^)nn) is given in Section VIA and 
shown in the inset of Fig. |[ This allows us to illustrate 
the SCHA in its simplest form. For this 2D Coulomb 
system we can approximate the equipartition result as 
(w^)nn ~ T/fit (i.e. we ignore compression modes, and 
assume no dispersion). The two equations relating (M^)nn 
and /Lit can then be solved graphically for any tempera- 
ture as in Fig. ^ and the maximum value of (u^)nn for a 
stable self-consistent solution is easily found. In this case 
we find the large value of (u^)nn/a^ ~ (0.34)^ at T„. 




FIG. 2. Illus tratio n of the SCHA for a 2D Coulomb sys- 
tem (see Section VIA). The inset shows the shear modulus fit 
of the 2D Coulomb lattice as a function of nearest-neighbor 
displacements (M^)nn, as given in Eq. (B9|). We show graphical 
solutions to the equipartition result, {u )nn ~ T/fit for three 
different temperatures. At high temperatures there are no 
self-consistent solutions. Note the non-perturbative form of 
fit at small fluctuations: The curve cannot be approximated 
by a power series in (u^)nn. In fact, further investigation of 
the diagrammatics shows that this result severely underesti- 
mates the thermal softening of the lattice. 



The following argument indicates that this instability 
point in the SCHA is a rigorous upper bound on the true 
stability limit of the crystal phase. If there is no solution 
to the SCHA then by definition there is no stationary 
value of the RHS of Eq. (|l^), which implies that this 
RHS is an unbounded function of trial parameters. Be- 
cause of the inequality, there can be no finite free energy 
F, and no equilibrium phase corresponding to the exact 
Hamiltonian H[uk]. Therefore the only possible equilib- 
rium phase is one without a well-defined mapping from 
particle positions to reference lattice points. In this case 
the Fourier transforms Uk are not the relevant variables 
of the system, and the liquid phase is still allowed. 

Unfortunately, it is known that the SCHA is not al- 
ways a reliable approximation, and we referred in the 
introduction to an example where it drastically underes- 
timates the thermal softening of the latticed To under- 



stand why this is so we must consider which diagrams the 
SCHA represents in a perturbation expansion around the 
harmonic limit, and this is done in Section IV B. 



IV. PERTURBATION THEORY BEYOND THE 
HARMONIC APPROXIMATION 



We now consider the standard perturbation expansion 
about the HAE3 The exact partition function is. 



Z = 



(17) 



where _ff'[uk] = i?[uk] — i?2[uk]- As usual we expand the 
exponential in 



Z=Zn 



(18) 



T 



With the same expansion, the propagator is, 
•' k' 

= «"^k>o (19) 

- (_l)n/ ^^^^ ^H'[V.^,] 

2^ T,I 



The brackets 



/c,Go 



denote a cumulant average in 



the HA, that is, using Wick's theorem the u fields must 
be paired up to give harmonic Green's functions. Go, and 
all terms with disconnected pairings are ignored, as these 
cancel exactly with the expansion of the partition func- 
tion Z in the denominator. To establish the terms in this 
series at each order in T, we need to separate each order 
in u of the Hamiltonian, 



-ff'[Uk] = i/3[Uk] + -ff4[Uk] + i/5[Uk] 



(20) 



where iJ,„[uk] contains all terms of order u™, and will 
lead to a vertex in the diagrammatic expansion with m 
legs. Exphcitly, 



"i! Jbz (27r) 



(ki 



(21) 



The crystal symmetry of the ground state lattice implies 
that the interaction tensor conserves pseudomomentum. 



A 



Ai...A„ 



(ki,...k„,) = (i"^) 



.6ui 



(22) 



= $^-^"(ki,...,k„,_i)A'*(ki 
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where A''(k) = ^''(k + Q^). The important point 
is that in the HA the expectation of the propagator 
(MkM-k)o is proportional to T, so that the average value 
of Hm is of order . We can write the perturbation 
expansion as, 



Ui,U' ,,) = (Ui^U 



-k/0 



(23) 




c,Go 



and then take all possible pairings using Wick's theorem. 
We can represent this series in diagrams and the terms 
to order are shown in Fig. ^. The order of each dia- 
gram (in T) is counted by adding up the number of lines, 
and subtracting the number of vertices. Notice that to 
order T^, only H-^ and H/^ contribute, while to order 
we also have diagrams including and ffg- This is a 
general rule: the terms -ff2m and _ff2m-i first appear to 
order T"\ 



(a) 



(b) 




(c) 



(d) 



(e)- 



FIG. 3. The diagrammatic expansion of the propagator 
as described in the text, to order T'^. We have drawn the dia- 
grams in four groups. Group (a) represents the expansion to 
order T^, which only involves third and fourth order vertices. 
The remaining groups are all of order T"^. Group (b) shows 
the reducible diagrams that are trivial to resum (see text). 
Group (c) shows the irreducible diagrams to this order that 
can still be generated by a self-consistent resummation of the 
lowest order diagrams. Group (d) shows the diagrams that 
involve only third and fourth order vertices, yet are topolog- 
ically new to this order. Group (e) shows the first diagrams 
that include fifth and sixth order vertices. 

We now make a resummation that allows us to express 
the full propagator G in terms of the harmonic propaga- 
tor Go and a self energy E. Note that, because wc do 
not include disconnected diagrams, each Uk in ( p3| ) must 
join with a Wk' in H' when using Wick's theorem. This 
means we can rewrite (p^ in the form. 



,13 



:) = ('ik"!k>0 

+ ««-\c>0("k 



(24) 



'^^-k>0 



OO 

E 

n=l 



{-If 



5u' 



k^"k 




c,Go 



Inspection of Fig. js] shows that many possible terms will 
be trivial "copies" of lower order terms. These are the re- 
ducible diagrams, for instance group (b) in Fig. ^, which 
have the property that when we group all of the reducible 
diagrams corresponding to a given irreducible diagram, 
they form a resummable series. We can write this in the 
compact form 

G"^(k) = Gf(]^) + G^^^(k)E^i^=(k)G^^'3(k), (25) 



with the self energy defined as, 



S-^i^^(k) ^TL-^Y^ 



(-1)" 



n=l 



~T 



(26) 



The subscript i in the average tells us to only include 
irreducible contributions, that is, any term which can 
be separated by cutting just one leg should be ignored, 
as it is already included in the full Green's function on 
the RHS of (p5[). A further resummation is possible by 
evaluating the RHS of ( ^^ with respect to the renormal- 
ized, rather than the harmonic, Green's function. The 
functional form of the self energy will then depend self- 
consistently on G, 



(k,G)= 



i-iy 



T 



(27) 



s,G 



where the subscript s tells us only to include "skeletal" 
diagrams in the average. We define a skeletal diagram as 
one that has not already been included in the resumma- 
tions of lower order diagrams by using the full Green's 
functions. For instance the diagrams of group (c) in Fig. ^ 
must be ignored, as they are included when we evalu- 
ate the self-energy contribution from group (a) with the 
renormalized propagator. After these resummations we 
have a Dyson-type equation, which must be solved self- 
consistently for G, 



(G-i)"''(k) = (Go')"''(k) - S"'3(k,G). 



(28) 



Note that this inverse propagator is exactly the response 
function <i>"^ (k) of the fluctuating crystal for an external 
force. This is demonstrated in Appendix 



A. The (3+4)-SCA 

Due to the difflculties of calculating all possible dia- 
grams in the perturbation expansion, one looks for ap- 
proximate resummations. A simple example is to only 
take the lowest order skeletal diagrams (group (a) in 
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Fig. W in the self energy, which is known as the (3+4)- 
SCAfl This is shown in Fig. ^, and it corresponds to 
solving the Dyson equation with the self energy, 



■^3+4 



(k,G) 




(29) 



Note that this method gives the exact propagator to or- 
der T^, although it only includes six of the thirteen dia- 
grams at order T^. The two skeletal diagrams are typi- 
cally of opposite sign and, importantly, for all the cases 
of 2D power law interactions looked at by Lozovik et 
al. the magnitude of the positive H^^ term is slightly 
greater than the negative term. This means the lat- 
tice first softens with temperature, but any approxima- 
tion that does not include the H^^ term will incorrectly 
find a hardening of the lattice with increasing tempera- 
ture. More details of the (3-|-4)-SCA and calculations for 
the example of logarithmic interactions in 2D are given 
in Appendix 







identical diagrammatic series appears for^a-form of SCHA 
in the context of the sine-Gordon modeled. In this case, 
the cosine potential only has terms even in the field, and 
so the SCHA has a better chance of being an accurate ap- 
proximation.) For our problem the SCHA only includes 
one of the two diagrams that contribute to the propaga- 
tor at order T^, and only three of the thirteen diagrams 
at the next order (compare to Fig. ^). The usefulness of 
the SCHA is that it is straightforward to calculate aver- 
ages of the true Hamiltonian over a Gaussian distribu- 
tion. This avoids the difficult task of actually evaluating 
high order vertices in a diagrammatic expansion. (E.g. 
for a pairwise interaction one must evaluate a sum over 
m derivatives of the potential to find the m-th vertex. 
This is done for m — 3 and to = 4 in Appendix p|.) 



FIG. 5. The diagrammatic representation of the SCHA. 
The thick lines are the renormalized propagators. The filled 
triangle represents the self-energy in the SCHA. The symme- 
try requirements of this approximation mean that only even 
vertices are included. 



FIG. 4. The diagrammatic representation of the 

(3-|-4)-SCA. The thick lines are the renormalized propaga- 
tors. This is just the Dyson equation for the two lowest order 
diagrams, the "Hartree" diagram coming from H4 and the 
"flying-saucer" diagram coming from (Hs)^. 



B. Diagrammatics of the SCHA 

It is insightful to consider which diagrams the SCHA 
(see Section [H) corresponds to. Writing H = H2 + H' 
we see that Eq. dl2) becomes 



(30) 



Comparing this to (28) and equating to a Green's 
function Gt~^ we sec that the SCHA corresponds to solv- 
ing the Dyson equation with the self energy given by 



■^SCHA 



(k,Gt) 



-L" 



(31) 



which is equivalent to the n — 1 part of the sum in (|27 
(for n = 1 there are no non-skeletal diagrams, so that the 
averages are equivalent). Therefore the SCHA is formally 
identical to the diagrammatic equation in Fig. Note 
that only even order vertices are present, as the aver- 
age of any odd combination over the elastic Hamiltonian 
is zero. This is an important feature of the SCHA. (An 



V. A NEW SELF-CONSISTENT METHOD 

We want to develop a method that combines the sim- 
plicity and elegance of the SCHA with the accuracy at 
lowest order of the (3+4)-SCA. We note that the SCHA 
is equivalent to evaluating just the n = 1 terms in the 
self energy (|6|) , while it is the n = 2 part that includes 
the "flying-saucer" diagram that is so important in the 
(3-|-4)-SCA. The contribution to the self energy from 
n = 2 is, 



{H'Y 



(32) 



s,G 



To write this as an average over the full Hamiltonian 
with a given distribution of displacements, we must sub- 
tract the diagrams that were included in resummations 
of the n = 1 diagrams. Allowing for this, and separating 
H' = H — H2 leads to the form. 



s;:L2^(k,Gt) = — 



1/ 

T 



5H SH 



(33) 



Note that terms explicitly including H2 have cancelled 
with each other. This corresponds to calculating, within 
an SCHA type method, the effective elastic matrix. 
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(34) 



We call this equation the "two-vertex-SCHA" . The ex- 
tra terms look rather like a fluctuation contribution to 
the elasticity, but we should remember that it is simply 
the two-vertex part of the Dyson series, and there will 
in principle be further corrections with higher numbers 
of vertices. The diagrammatic representation of the two- 
vertex-SCHA is shown in Fig. ^ This scheme now incor- 
porates the diagrams from the SCHA and the (3-|-4)-SCA 
along with more terms not included in these last two ap- 
proaches. For instance, it includes nine of the thirteen 
diagrams for the propagator at order T^, see Fig. ||. It 
also has the feature of the SCHA that one only needs to 
calculate averages of derivatives of the full Hamiltonian, 
and not any diagrams explicitly. 




FIG. 6. The two-vertex-SCHA described in diagrams. 
This approximation includes nine of the thirteen diagrams 
to order T"' in Fig. ^. This is equivalent to the self-consistent 
solution of Eq. (|34[), so that the diagrams do not need to be 
explicitly calculated. 



VI. PHYSICAL REALIZATIONS 

A. The 2D Coulomb system 

The 2D Coulomb system, also called the one- 
component plasma, consists of particles interacting with 
a pairwise, repulsive logarithmic potential, 



V{R) = -uln(i?/C). 



(35) 



It has been studied in the pastp^artly because of some 
special mathematical propertiesEZl (the system is exactly 
soluble at one particular temperatuce pi= v/2 in the liq- 
uid phase). Numerical simulationaHaL^ have found evi- 
dence for a weak first order transition at ~ 0.007?;. 
Physically, this system is of interest as it is applicable to 
a thin superconducting film ysdth field-induced vortices 
which interact logarithmicallyEil up to a screening length 
A = A^/d where d is the thickness of the film and A is the 



bulk penetration depth, typically of order 1000 A. Such 
films were predicted to undergo a 2D melting transition 
of the vortex Vjttice well below the zero-field transition 
temperatureEj'tJ (A similar scenario takes place in lay- 
ered superconductors at high fields, see Ref. [lO| .) 

The ground state of the 2D Coulomb model is a 
triangular lattice, with nearest-neighbor spacing a = 
(2/nV3)^/^. This lattice has a shear modulus of 



Mo 



lim 



(36) 



(this result was known over thirty years agccj) and a 
diverging compression modulus at large wavelengths. 



Ao(fc) = $§"(fcx)/fcx' ~ 2nvn'/e. 



(37) 



In Refs. ^ and ^ it was presumed that the melting tran- 
sition would be of the continuous Kosterlita-Thouless- 
Halperin-Nelson- Young (KTHNY) typeyfl with a 
dislocation-unbinding mechanism, which is predicted to 
occur at 



Tkthny — A* 



47r 



A{T) 



16%/37r 



O.Oll^(T)?; (38) 



with A < 1 the renormalization constant for the shear 
modulus at finite temperature, /i — fi()A{T). How- 
ever, this iS|_aQi_consistent with the first-order result of 
simulations .Ea^Lj One possible reason for this discrepancy 
is that regular fluctuations in an anharmonic potential, 
such as we considered in the last section, can destabi- 
lize the lattice before singular fluctuations such as free 
dislocations are generated. Therefore in this section we 
consider the anharmonic softening and instability of the 
2D Coulomb lattice, and compare to the predictions of 
the KTHNY theory. We calculate ^ within the SCHA. 
We also find the first order correction in temperature to 
the shear modulus using perturbation theory, by calculat- 
ing the diagrams of group (a) in Fig. |3[ This shows that 
the SCHA greatly underestimates the thermal softening. 
We then present our results for ^ from the (3-f4)-SCA 
and the two-vertex SCHA and compare to the KTHNY 
universal value of /i = 47rTKTHNY / ■ 

The SCHA for the 2D Coulomb system was briefiy 
looked at in Ref. |2^. In fact, the SCHA with "indepen- 
dent fiuctuations" has some nice features for this system 
due to the Laplacian of the interaction being zero. In 
Appendix ^ we derive the non-perturbative form of the 
elastic moduh as a function of a^o — i(|u(R^) — u(0)|^), 
which leads to the result for the shear modulus, 



MO X! ( 1 

« /xo + 6^0 ( 1 



R 



2a 



(39) 

aV2(«')n„^ (40) 



This formula is shown in the inset of Fig. ^ (the nearest- 
neighbor approximation here is extremely accurate). The 
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non-perturbative feature (the extreme flatness at small 
(w^)nn) will be lost when we include the correlations in 
fluctuations between different particles. In fact, we have 
found that the SCHA predicts an initial stiffening of lat- 
tice with small fluctuations, which we understand to be 
an incorrect result due to the neglect of diagrams with 
odd vertices. 

Inspection of Eq. ( ^ ) shows the correct temperature 
expansion of the change in shear modulus to be, 



a = tin — lim 



(ky, G)/ky 



(41) 



We can find the first correction to the shear modulus 
to order T by evaluating the bare self-energy contribu- 
tions in the diagrams of Fig. ^ Formulas for the (3-1-4) 
self-energies for a system with pairwise interactions are 
given in Appendix O We have numerically evaluated 
these formulas for the 2D Coulomb system, and our re- 
sults for the transverse Yj^ = I]^^(fca;) and longitudinal 
= Ys'^[kx) parts taking G — Gq are shown in Fig. |^. 
The large wavelength limits are found to go as, 

j:J{k) w {UT/a^)k^ (fc) « -(10r/a2)fc2 (42) 



I]f'(fc) « 98T/a^ El'(fc) « {lQT/a^)k'^ 



(43) 



100 R 




FIG. 7. Contributions to the transverse and longitudinal 
self e nerg y from the skeletal diagrams of Fig. |^, as given by 
Eqs. (D3) and (D4), for logarithmic pairwise interactions in 
2D. 



From Eqs. (41 ) and ( f42| ) the shear modulus at low tem- 
peratures will be given by 



/i w ^0 1 - 28 



T 



(44) 



so that within first-order perturbation theory, /i is about 
80% of its zero temperature value at the numerically ob- 
served transition at ~ 0.007w. The large wavelength 



limit of the longitudinal part of the self energy in (|42 
is a constant, so that the compressional modes are also 
softened to first order in temperature. 



T 
12- 

V 



(45) 



This thermal softening of the long-wavelength compres- 
sion energy (a special feature of the 2D Coulomb system) 
might seem to contradict the physically intuitive expec- 
tation that the compressibility is fixed by the bare charge 
of the system. A full understanding of this result will be 
discussed in a future paper, but it is due to the non- 
linear relation between changes in the density and the 
longitudinal displacements, Sn = nV • u + 

We can go beyond the lowest order in T using the 
(3-1-4) self energy by iterating the equation of Fig. ^ 
until a stable solution is found within this (3-f4)-SCA. 
We have also found numerical solutions within the two- 
vertex-SCHA of Section 0. Details for calculating the 
effective elastic moduli of Eq. (|34| ) for pairwise interac- 
tions are given in Appendix K[ ^As we have shown in 
Section 0, this approximation gives identical results to 
lowest order in T as the (3-|-4)-SCA. However, we find 
important deviations at larger temperatures. In Fig. ^ 
the results for fj, as a function of temperature from both 
methods are shown. The shear modulus was recently 
measured in—the same system with Langevin-dynamics 
simulations ,E3 and we find good agreement between the 
measured temperature dependence (the circles in Fig. |^) 
and the two-vertex-SCHA calculation. Notice that the 
instability point from the (3-|-4)-SCA is below the nu- 
merically measured melting temperature of T,„ ~ O.OOTt;, 
while the upper limit to solutions of the two-vertex- 
SCHA is slightly above this temperature. 

In Fig. H we also compare to the predictions of the 
KTHNY theory of 2D melting. As described above, this 
theory predicts a universal value of the shear modulus at 
the dislocation-unbinding transition, — At:T /a? , shown 
as the dot-dashed line in Fig. ^. The fact that the two- 
vertex-SCHA result crosses this value tells us that the 
results at temperatures above this value are unphysi- 
cal, as the crystal is unstable to free dislocations (if the 
(3-|-4)-SCA were correct, it would practically rule out 
the KTHNY mechanism for this system). Therefore the 
crossing of the universal jump and the two-vertex-SCHA 
result is an upper bound to a melting transition and we 
find Tu — 0.0072w (the effect of thermally generated dis- 
location pairs will push /i, and the melting temperature, 
down further). The fact that at melting the shear mod- 
ulus is so close to the KTHNY universal value tells us 
that dislocation-pairs (not included in our analysis) as 
well as anharmonicities both play an important role in 
the physics of the melting transition. Whether the weak 
first-order nature seen in simulations is correct, or a nu- 
merical artifact, remains unclear. 
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FIG. 8. The shear modulus of the 2D Coulomb crystal as a 
function of temperature as calculated to first-order in Eq. ( jji] ) 
(short dashed line), and within the (3+4)-SCA (long dashed 
line) and the two-vertex-SCHA (full line). At higher orders 
in T we see that the two-vertex-SCHA predicts less softening 
than the (3-|-4)-SCA, and a correspondingly higher instabil- 
ity temperature, T^. The circles are the measured values of 
/i from numerical simulations taken from Ref. 36. They agree 
well with the two-vertex-SCHA curve. The dot-dashed line 
follows the universal jump prediction of the KTHNY transi- 
tion, [i = A-kT/o?. Below this line the lattice is unstable to 
free dislocations. 



B. Classical electrons in two dimensions 

It is experimentally feasible to confine electrons to a 
flat surface in various systems, such that they are free to 
move in only two dimensions. In particular, for electrons 
trapped on the surface of liquid helium by an applied 
electric potential, a classical melting transition from a 20, 
"Wigner crystal" to a liquid state has been observed.E3 
The elasticHproperties of the 2D Wigner crjgtal are well 
studiedJ23'L3 and Platzman and FukuyamatI considered 
the effect of thermal fluctuations using the SCHA, report- 
ing a dramatic stiffening of the lattice with temperature^ 
Using the perturbation expansion of Section FisherO 
showed this to be incorrect, as the first order term in T 
gives a downward correction to /i. Lozovik et al. have ap- 
plied their (3-|-4)-SCA to this problem,EI and stated that 
an anharmonic instability is reached below the KT HNY 
temperature. However, the discrepancy in Section VIA 



for logarithmic interactions between the (3-|-4)-SCA and 
our two-vertex-SCHA motivates us to reconsider the 2D 
Wigner crystal. 

Electrons interact with a repulsive pairwise potential. 



V{R)^ 



R' 



(46) 



where e is the electron's charge. The ground state is a tri- 
angular lattice, spacing a, and at finite temperatures all 



properties depend on the dimensionless ratio T/(e^a ^). 
In the literature the parameter F = ^Jlin'^e^ jT is more 
often used. When calculating the elastic matrix, care 
must be taken to split the interaction to long- and short- 
range parte, and treat them separately as in the Ewald 
techniqueJ23 Then the zero-temperature shear modulus 
is found to be. 



Mo 



0.2450645 e^n^/^ 



(47) 



and the compression modulus diverges at small k (as for 
the Coulomb crystal but with a different power) as 



Ao(fc) 



27re^n^ 



(48) 



Repeating the perturbation theory of Section IV we find 
the first-order correction to the shear modulus to be. 



T 



(49) 



in agreement with Ref. |g. This correction again comes 
from a negative contribution from S|"(G'o) and a positive, 
but smaller in magnitude, contribution from ^^(Go). We 
have also confirmed that the SCHA gives an incorrect 
stiffening of with increasing temperature. 
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FIG. 9. The shear modulus of the 2D Wigner crystal as a 
function of temperature. We show the first-order correction of 
Ref. 8 (short dashed line), and our results for the (3-|-4)-SCA 
(long dashed line) and the two-vertex-SCHA (full line). The 
results are qualitatively the same as for the Coulomb crys- 
tal in Fig. H with the two-vertex-SCHA result crossing the 
KTHNY universal jump line (dot-dashed line). 

In Fig. ^ we plot our results for the shear modulus /i(T) 
from the (3-h4)-SCA and the two-vertex-SCHA. The re- 
sulting picture is very similar to that for logarithmically 
interacting particles in Fig. ^. Our instability point for 
the (3-h4)-SCA of r3+4 « 155 is close to the quoted 
value from Lozovik et al. of F^"'"'' = 140. The two- 
vertex-SCHA is again more stable, and together with the 
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KTHNY condition on the minimum allowed ^(T) gives 
~ 135. This is close to ptii|e-|experimental melting 
temperature of T„i = 127 ± 5j23'c2l valid over a range of 
experimental densities. Also, numerical simulations have 
found a melting transition in the range from F^™ « 90 
(Ref 
(Ref. |44|) 



A 



L|) to r^™ « 125 (Ref. ^ and |43D to T^™ « 159 
(It is not clear why there is such disagreement 
between the numerical results. The rathetJiigh melting 
temperature from the earliest simulatiorc2l is probably 
an out of equilibrium artifact; the system was started in 
the crystal state and heated up.) 



C. Quantum melting of 2D bosons and the Bose 
model of vortex-lattice melting 

In the previous sections we have ignored the effect of 
quantum fluctuations. The uncertainty principle can pre- 
vent the formation of a lattice phase if the mass of the 
particles is sufficiently low. In this section we will con- 
sider the effect of quantum fluctuations on a 2D Coulomb 
lattice of bosons at zero temperature. In the path- 
integral representation, the action of the boson system 
in imaginary time is, 



S 



h/T 



dr 



E 



dK,, 



dr 



R,, 



(50) 



with V{R) given in Eq. (p35|). The properties of the sys- 
tem are determined by the statistical sum, or partition 
function, Z = TrR^(^) [e"'^/'*] . At zero temperature all 
properties depend only on the unitless de Boer param- 

1 /2 

eter, A — (fi^ / a^vra) , which measures the ratio of 
zero-point kinetic energy to the interaction energy. 

This model maps directly to a 3D classical system 
at finite temperature, where the world lines of the 
bosons correspond to fluctuating elastic strings. In fact, 
the Coulomb interaction between the strings provides a 
model of the interactions between vortex lines in a type- 
II superccmductor (this fruitful analogy was discovered 
by Nelsonc3). Thus the path- integral representation of 
2D bosons has been studied with, the motivation of de- 
scribing vortex-lattice meltingi3'L3 which oceurs at finite 
magnetic fields in high- Tc superconductors.lly To be spe- 
cific, the 3D vortex system is obtained by replacing h 
with T , the time r with the distance along the field di- 
rection z and the mass m with the line elasticity (In 
fact, this model only approximates the 3D interactions 
of the vortex system. A more correct mapping would be 
from a boson system with interactions non-local in time.) 
The strength of the interaction per unit length is deter- 
mined by the penetration depth A of the superconductor 
via the relation v 2eo with Eq = ((/)o/47rA)^ where the 
quantum of flux is c/jq — hc/2e. In these units the de Boer 
parameter is 



(51) 

In the numerical simulationsllBS a melting transition oc- 
curs at A « 0.062, above which the system is in a line- 
liquid phase. Although it has proved difficult to mea- 
sure the equilibrium values of the elastic moduli in the 
same simulations, the value of (u^) was measured, and 
shows deviations from a linear dependence on tempera- 
ture before melting occurs.CJ In this section we use our 
self-consistent treatment to find the effect of anharmonic- 
ity and the instability point in this system and compare 
our results to the measured (u^) in the simulations. We 
also calculate the quantitative effect of fluctuations on 
the shear modulus of the vortex crystal. 

In the following we use the notation for the classical 
vortex system, but all results are easily transferred to the 
2D quantum system at zero temperature. The general- 
izations to fluctuating lines of the results of Sections 2-5 
are straightforward, where the displacement variables are 
2D vectors in 3D k-space. 



1 



Uk 



d^y^ u^(z)e" 



-i(kj^.Rj,+/c2z) 



(52) 



The elastic matrix and the self energy are now functions 
of kz as well as the 2D vector kj^. Brief details of the 
(3-t-4)-SCA and the two-vertex-SCHA are given in Ap- 
pendix ^. 

At large wavelengths the transverse elastic matrix is 
given by lims-^o $"^(k) — CQ^k^ + 044^^ and the longitu- 
dinal modes by limfc^o $^(k) = ciik']_ -1-044^^ where cn, 
C66 and C44 are the Voigt notation for the compression, 
shear and tilt moduli respectively. We can again find the 
first correction to the shear modulus by evaluating the 
self energy to lowest order in T. We find that 'S'^{k±) w 

3.2(T/a2a^)fc2^ and T.J{kj_) w ~2.0{T/a^az)kl, where 

1 

~ a (£;/eo) ^ . Therefore to lowest order in the de Boer 
parameter, 



C66 ~ (1 - 5.5A) , 



(53) 



where the shear modulus without fluctuations is Cgg = 
neo/4. The compression modulus at large wavelengths 
is not renormalized. Note that the "Hartree" self-energy 
S4(k) is independent of the z-component of the wave- 
vector (see Appendix^), but the contribution S3(k) does 
depend on k^. For example, for k^ at the Brillouin zone 
we find, Y.f{kl = kBZ,kz) w (T/aV)(8.9 - 2.8alkl), 
which gives a fiuctuation induced stiffening of the tilt 
modulus. 



C44(fc5 



k 



BZ) 



El 

«2 



l-t-3.9A). 



(54) 



We have made detailed calculations for the two-vertex- 
SCHA for this model, and the results for the shear mod- 
ulus as a function of A are shown in Fig. It turns 
out that the effective elastic matrix (see Appendix ^) 
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depends on the full correlations of displacement fluctu- 
ations of a given vortex line in the z-direction, and we 
need to numerically find a self-consistent solution for the 
function, 



= T 



u"(z)u"(0)) 

, 2lT 



(55) 



BZ 



=(i>r')""(k). 



In practise we discretize along the z-direction, making 
sure that (m^(z)) is smoothly varying. 
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FIG. 10. Shear modulus for 2D bosons with quantum 
fluctuations, or for the Bose model of the vortex lattice at 
flnite temperature. The long-dashed line is the first-order 
result (pS). The short-dashed line is the result from the 
(3-|-4)-SCA, while the full line is the shear modulus calcu- 
lated in the two-vertex-SCHA. Note that a stable solution is 
found for much smaller values of the shear modulus than for 
a two-dimensional system at finite temperature (see Fig. 

We again see that the stability limit from the two- 
vertex-SCHA is higher than that from the (3-|-4)-SCA. 
Surprisingly, the crystal maintains a self-consistent solu- 
tion for a shear modulus which is only 20% of the zero- 
temperature value. This contrasts with our results for a 
2D crystal, and is mainly due to the stiffness of the elas- 
tic lines: the average displacements within the equipar- 
tition theorem only gmw with the inverse square-root 
of the shear modulus,Ej [v?) cx T/^C66C44. The result 
is qualitatively similar to the measured shear modulus 
in simulations of coupled vortex pancakes in a layered 
superconductor in Ref. In their somewhat different 
model of a vortex lattice, a reduction in the shear modu- 
lus at melting of about 30% of the zero temperature value 
was found. Also, we find that the tilt modulus stiffens 
slightly with the fluctuations, allowing for a still softer 
shear modulus before instability is reached. 
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FIG. 11. Comparison of the single-particle fluctuations as 
measured in Monte Carlo simulations of Ref. 46 and as cal- 
culated in this paper. The circles are data from Ref. 46 and 
the full line is from the two-vertex SCHA, with a stability 
limit close to the numerical melting. The long dashed line is 
the harmonic result. The (3-f4)-SCA (short-dashed line) be- 
comes unstable well below melting, while the instability point 
of the SCHA (not shown) overestimates the melting point by 
a factor two. 

In Fig. nj we compare the mean fluctuation per parti- 
cle {u'^(Oyf~= {2n)~^ J (fk{u^ut^ with the values mea- 
sured in simulations (see Fig. 8 of Ref. ^6|). The sim- 
ulation results come from fitting a Debye- Waller factpt 
to the measured widths of the crystal structure factor.E^I 
The main point is that the instability point predicted 
by the two-vertex-SCHA, A^ = 0.061 is very close the 
numerically measured melting point A™ = 0.062. The 
fact that we find A„ slightly below melting means that 
the two-vertex-SCHA still slightly overestimates the level 
of softening, and correspondingly there is an increasing 
discrepancy with the numerical results for {v?) close to 
the melting temperature. We should of course remem- 
ber that the two-vertex-SCHA is only a limited resum- 
mation of the perturbation theory. However, our new 
method clearly gives much superior results compared to 
the (3+4)-SCA, also shown in the figures, and to the 
SCHA, which we found to give very little change to the 
elastic moduli below the numerically found melting and 
an instability bound over twice the value of Km- 

The thermal softening of the shear modulus of the vor- 
tex lattice is rarely explicitly calculated in the literature. 
The measurements from simulations in Ref. |3^ agree with 
Fig. |l^ in that this softening is really significant. It may 
help to explain the peak effect in the critical current seea 
close to the vortex-lattice melting transition in YBC0.E3 
The weak collective pinning theory predicts that, for a 
non-dispersive tilt modulus, the critical current depends 
on the elastic moduli as ici-y i'^/c^qCj^a where 7 is a 
pinning strength parameter.cZI However, the effect of fi- 
nite temperature is not only to reduce the shear modu- 
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lus, but also to smooth the effective pinning landscape 
over (u^), thus reducing the pinning force. To see if our 
results lead to a critical-current peak as a function of 
temperature requires a careful analysis of the k depen- 
dent elastic moduli (the critical pinning volume depends 
on long- wavelength distortions but (u^) comes from in- 
tegrating over all modes), and we leave this for a future 
project. 

VII. CONCLUSIONS 

In this paper we have introduced a new method, the 
two-vertex-SCHA, to treat the anharmonicity of simple 
lattices in the presence of thermal, or quantum, fluctu- 
ations. The method is based on the diagrammatic ex- 
pansion about the HA, but takes the form of an exten- 
sion to the SCHA calculation of effective elastic moduli. 
By considering two long-range interacting systems at fi- 
nite temperature, and a system of zero-temperature 2D 
bosons, we have seen that the two-vertex-SCHA is a sig- 
nificant improvement over previously used self-consistent 
techniques, the SCHA and the (3-|-4)-SCA, which respec- 
tively under- and overestimated the fluctuation-induced 
softening. It is the comparison to simulation results that 
gives us the most confidence in the reliability of the two- 
vertex-SCHA as a tool for calculating the softening and 
instability of such simple crystals due to thermal or quan- 
tum fluctuations. One should therefore not trust the re- 
sults of previous, treatments of the anharmonic softening 
of 2D lattices,LlEl but use instead the two-vertex-SCHA. 

In addition to establishing this new method, we have 
found important results for the vortex-lattice melting 
problem, using the Bose model for the finite-temperature 



vortex system. We have found a much larger than ex- 
pected thermal softening of vortex lattice with the shear 
modulus reducing to 20% of its zero temperature value at 
the melting transition, which may have important conse- 
quences for the pinning properties of the vortex lattice. 
It should be noted, however, that the real vortex system 
does not only have interactions within planes of equal 
height, as in the Bose model. The full London model 
of vortices correspondiLto a 2D Bose system with time- 
nonlocal interactionsO An investigation to extend the 
results of tiie Bose model to the full time-nonlocal case is 
underway.EJ Another facet of the vortex system in high- 
Tc superconductors is the layered nature of the cuprates. 
For very weak inter-layer coupling, the system can be de- 
scribed by magnetically coupled "pancake" vortices. In a 
recent paper the two-vertex-SCHA has been used tOj^d 
the instability line of such a pancake-vortex lattice.o It 
lies just above the melting line as calculated with a free 
energy comparison using results from 2D simulations for 
the pancake-vortex liquid phase. Finally, it would also 
be interesting to investigate the stability of lattices in 
the presence of quenched disorder. This is also relevant 
for the vortex lattice, where pinning disorder exists, and 
a disorder-induced melting transition occurs as the vor- 
tex density, or magnetic field, is increased.!^ 
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APPENDIX A: NEAREST-NEIGHBOR FLUCTUATIONS IN TWO DIMENSIONS 

In the SCHA defined by Eq. (|lj), we have to take averages over the second derivative of the potential between two 
particles as they fluctuate within the crystal (similar averages are needed in the two-vertex-SCHA, see Appendix ^). 
Within the Gaussian approximation these averages only depend on the mean square, a^o = 5(|u(R^) — u(0)|^), which 
in two dimensions is given by 



T 



BZ 



(1 



„ikR„ 



1 



1 



T 

ev 



In [R^/a) 



(Al) 



{ev is an energy-density scale set by the interaction V). Note that the fluctuations of a single particle, (u^) — 
lim/j^^oo(cr^o) diverges in 2D. We show here that the logarithmic divergence at large distances is not important for 
the elastic moduli which have dominant contributions from fluctuations at small distances. Consider the difference 
between the renormalized and the zero temperature elastic matrix 



$f (k)-$^"(k) 



n ^(cosk ■ 



R,, 



1) 



fq'^V{q)e 



j9 CTjio 



(A2) 



Now the Fourier transform only has weight at values of g < so as long as a^o << Rj^ we can write e 2? '^^o 

1 - iq^o-^o and, 
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^f{k) - <^f{k) « n^(cosk • - l)-a^oy^dadfiViR^) (A3) 

Now the dependence of cr^io is logarithmic and weak compared to the short-range variation in (cosk • — 
1)'S/^ dad/^V (Rf^) , so short distances will dominate and we are justified to approximate fluctuations by the nearest- 
neighbor result (Tpo ~ {u'^)nn- Similar arguments apply for the sums that appear in the two-vertex SCHA. (For 
logarithmic interactions, W'^V{R) — 0, and the perturbation theory is not helpful. However, for this case the non- 
perturbative result ( ^9| ) derived in Appendix |^ is clearly dominated by the nearest neighbors.) 



APPENDIX B: SCHA CALCULATION OF THE FLUCTUATION-RENORMALIZED ELASTICITY FOR 

LOGARITHMIC INTERACTIONS IN TWO DIMENSIONS 

In this appendix we show how to calculate the elastic matrix of a 2D lattice with logarithmic interactions within the 
SCHA as given by (14) and demonstrate the non-perturbative dependence on the size of fluctuations (u^). We need 



the thermal average of the second derivatives of the interaction potential between each pair of particles. Assuming 
the displacement fluctuations to be Gaussian distributed, this is, 

52F(R^ + u,,-uo)\ f d^rd^V{R^ + r) _J' nnji ^ ^\ 



dR'^dRP / J Z dr'^drP ^ ' dR^'dRi 
where the derivatives on the RHS do not act on the distribution function, 

V{v)=e^v[^^\r^K^r^rPy (B2) 



= (K(R^) - ^"(0)]K(R^) - ^'^(0)]), (B3) 



The width of the fluctuations is given by the matrix 

afi 

and Z = 27r[det(cr;^^)]5. 

We take the isotropic fluctuation approximation, a"^ = Sapo^Q, which we find to be good from numerical evaluation 
of (^ in the HA. In this case, because the potential is circular symmetric, the average will be independent of the 
angle of the ground state position R^, so 



27rJ Z 



inR,)) =7;z I —V{r + -R^)e -.o (B4) 



f d^r f " d 

^ / d'pWlVir + p)e 

-n- J p<R 



27r . 

1 f d^r /"^^ dR 
2^. ~Z 



where we have made use of the divergence theorem J d9^ (d/dR) = R ^ § dlh-V = R ^ jp^p d^P^% We are inter- 
ested in the case of logarithmic interactions, V{R) = — ?;ln(i?/^), and we use the fact that V^[ln(i?/^)] = 27r(5^(R) to 
give 

^ e-«V2-.o _ I ^ -v\n{R^/0 - V / ^g-^ /^a.o + ^onst, (B5) 
R L -I Jr.,, R 

where the integration constant depends on cr^o smd ^ but not on i?^. Notice that the correction term to the zero- 
temperature potential is non-perturbative in the fluctuation tr^o- The physics of this result is clear: because the 
2D-Laplacian of the interaction potential is zero, Gauss' theorem applies and the circular symmetric potential gives 
the same result as if all of the "mass" was at a point in the center of the circle (the ground state position R^^). The 
non-perturbative correction appears from the exponentially small tail of the distribution which is outside of the radius 
i?^. This "shell" does not contribute to the average interaction, so it's weight must be subtracted from the bare 
interaction. We can easily differentiate the correction term with respect to the lower limit of the integral, to find 
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- uo) \ _ d^V{R, 



dR'^dRP 



Ru^ 



R 



5a[: 
Rn 



(B6) 



This expression should be inserted into (jlj) to give the renormahzed elastic matrix. It is then straightforward to 
derive the softening of the shear modulus within the SCHA as given in Eq. (H) and Fig. | 



APPENDIX C: EQUIVALENCE OF ELASTIC RESPONSE WITH INVERSE PROPAGATOR 

We demonstrate in this appendix that the elastic response to an external force is given by the inverse propagator 
(or Green's function). In the presence of a force acting on each particle, the Hamiltonian is changed to, 

Hf[u,]^H[u^]-J2ff.-^,- (CI) 
The definition of elastic response within the thermally fluctuating region is that the force is given by Hooke's law, 

f^ = Y,^-P(^R^^n^)u^^ (C2) 

which is valid to lowest order in the mean displacements IJ^ — {u^)f- The displacements are calculated with 
straightforward statistical mechanics. 



If we have lattice-translational symmetry, the Fourier transform of this reads as = f^G"^(k), where the Green's 

function is defined by {u^u^'} f=a = TG'"''(k)5(k + k'), so that the elastic matrix in k-space is exactly the inverse 
Green's function, 

^'"''(k) = (G-i)"''(k). (C5) 

Note that this example of the classical fluctuation-dissipation theorem goes beyond the equipartition result, which is 
only valid for the harmonic approximation. The true elastic response of the fluctuating system is always given by the 
inverse of the propagator. 

If we construct the free energy corresponding to Hf, we will find the result, 



Ff = Ff^o ~\l 7f^(G-^)"^(k)L/,"C/f , + 0{f). 



(C6) 



This is for a system with forces specified. However, one is often interested in the case where a certain displacement 
is enforced (e.g. when boundary conditions are changed in a numerical simulation). We must then make a Legendre 
transformation to the relevant thermodynamic potential that depends on a fixed displacement, 

A(T,[U^])=F/+^/;{/; (C7) 

^Ff=o + U ^{G-'rf^ik)UQU^_^. + 0{U'). (C8) 
^ Jbz (^""j 

With this formulation, we see that the elastic matrix may be found as the second derivative of the thermodynamic 

potential with fixed displacements, $"'^(k) = d"^ A/ dU^dUty, . For example, the A; ^ shear modulus may be 

measured as the response to the change in the angle 9 determining the boundaries of a finite system as in Ref. ^ As 
the Hamiltonian will depend on the "displacement" 6*, we can differentiate the relevant free energy, A{T,9), to give. 



dO^ /g=o T\[de) ]_^^T\d6/,^, 



9=0 



(C9) 



The important point is that this shear modulus is exactly the same as the shear modulus found from the fc — > limit 



of the inverse of the transverse propagator ^ = hmfc^^o[(G ^Y^{ky)lky]. 
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APPENDIX D: THE (3+4)-SCA FOR PAIRWISE INTERACTIONS 



In this appendix we show how to exphcitly calculate the (3+4) diagrammatic expansion of Fig. ^, for the case of 
pairwise interactions. For a Hamiltonian of the form the third and fourth rank anharmonic tensors defined by 
Eqs. @ and (|2|) will be given by, 

$^i^^^3(ki,k2) =„y^ fe'ki.R,, ^g»k2.R^ _^g-*(ki+k2).R,A ^ 



J2 /3(Qm + ki) + /3(Qm + k2) + /3(Qm - ki - k2), (Dl) 



and 



$^i^=^-^^''(ki,k2,k3) = (1 - e'*'!-^") (1 - e'''^-^") (l 



/4(Qm) - /4(Qm - kl) - /4(Qm - k2) - /4(Q^ - kg) - ^(Q^ - ki - ks - kg) 
+ /4(Qm - kl - k2) + /4(Q^ - k2 - ka) + /4(Q^ - kg - ki), (D2) 

with /m(Q) = Q\i ■•■Qa„^(Q)- The above Poisson resummations to reciprocal-lattice sums allow for more con- 
venient evaluation, especially in the case of long-ranged interactions. The contribution to the self energy from the 
"flying-saucer" diagram is. 



S^'^(k,G) = — 



while the contribution from the first, "Hartree"-type diagram is 



? / ^G^^'*(ki)G^^^nk + ki)<i>^i^3^=(k,ki)$^^^*^«(k,ki), (D3) 
2 Jbz (Stt)" 



E^^^(k,G) = -LV^-^ ) ^G^^^^(ki)ci>^i^^^^^^(k,-k,ki). (D4) 



APPENDIX E: THE TWO-VERTEX-SCHA FOR PAIRWISE INTERACTIONS 

We now describe the details of the two-vertex-SCHA of Section ^ for a system with pairwise interactions. Splitting 
the efi^ective elastic matrix of Eq. (pj) term by term we have. 



$r''(k)-$f (k)+ci>f(k)-fci>^f(k), 
with $"'^(k) given by the SCHA formula (|l|), 



(El) 



$f (k) = - 



SH SH 



and, 



-1 



E 

a 



ik R,, 



ikR„ 



zk-Ro 



/ dV 


dV 








p-a) 



(E2) 



<(k) 



/5A2 



^-k"k 



(E3) 



where the averages are taken with respect to the effective elastic Hamiltonian. In the independent fluctuation approx- 
imation, we can write the entire elastic matrix as a function of ( u^) and T (note that in 2D, although (w^) diverges, 
we can take the nearest-neighbor value (w^)nn as the sum in (E2) is dominated by short distances). 
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^ ^ V- r a/3 _ a/3 _ a/3 a/3 

^2 W — J, [yQ,Q' yQ-k,Q' yQ,Q'+k + yQ-k,Q'+k 



Q,Q' 



n2_ r_A^V- r a/3 _ ^a,9 _ , a/3 

(27r)'^'^ r Q-q'-k,q' '*Q-q'+k,q' 

Q 



with 



and 



Finally the correction term is given by, 



-qi.q2(«^)/<i 



-qi.q2(ti')/d 



- 1 



$2f(k)-£^<i>f (k)$f(k). 



(E4) 
(E5) 

(E6) 
(E7) 



To complete the two-vertex-SCHA one must find a self-consistent solution numerically for {v?) by inserting the effective 
elastic matrix in the equipartition result (^. 

APPENDIX F: MODIFICATIONS FOR INTERACTING ELASTIC STRINGS IN THREE DIMENSIONS 

We now outline the generalization of the self-consistent methods to the problem of elastic strings directed along z 
and interacting within 2D planes of equal z. First, the harmonic energy is given by, 



Ho = 



1 



BZ 



(27r) 



27r 



^k"-k- 



(Fl) 



The generalization of Eq. (D3) is 



and of Eq. ( p4| ) is 



I]^^^(k,G) = -^ iez^ I ifG^^'^(k')*2V^'^'^(k^,-kx,k', 



(F3) 



The matrices with 2D subscripts are to be calculated for the interactions within a given 2D plane. 

For the two-vertex-SCHA we now have to be careful to include the z-correlations in the fluctuation terms, — 
(m"(z)w"(0)). We can again split the effective elastic matrix as in (El), with $i being the 2D SCHA elastic modulus 
plus the single string elasticity, and 



$f (k) = ^ E 

Q Q' 



dz e-'^'"' 



2T J (27r)2 



E 

Q 



oh"'^ (z) - h"'^ (z) - h"'^ (z) 

^"■Q-q'.q'V-^^ "Q-q'-k^qn^^ ""Q-q'+k.qH^ J 



(F4) 



dz e~*'= 



with 



and 



-iqi.q2(ti^(z)) _ Y 



where {u'^{z)) is defined in Eq. (p5|). The last, correction term is given by, 



<(k) = ^<i>f2D(kj.)<i>f2Z5(k±) / dz, 



— ik^z ,2 



{z)). 



(F5) 
(F6) 

(F7) 
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